package src.montecarlo;

import java.io.File;
import java.io.FileNotFoundException;
import java.io.FileOutputStream;
import java.io.IOException;
import java.util.ArrayList;
import java.util.List;
import java.util.Random;

public class MetropoliMonteCarlo {
	
	private int sampleSize;
	private double delta;
	
	Random gen;
	
	List<Double> results;
	
	public MetropoliMonteCarlo(int sampleSize,double delta){

		this.sampleSize = sampleSize;
		this.delta = delta;

		gen = new Random();	
		results = new ArrayList<Double>();
	}

	/*
	 * returns number of rejected trials.
	 * generates results on output file
	 */
	public int metropolisAlgorithm(String outputFile) throws IOException{
		
		//clear, just in case some previous results were stored
		this.results = new ArrayList<Double>();
		
		int numberOfSamples = 0;
		int rejected = 0;
				
		double x = 0;
		double previousProba = 0;

		//initial trial step: ideally should be close to the probability maximum
		do{
			//drift from the previous position
			x += (2*gen.nextDouble()-1)*delta; 
			
			double newProba = normal(x,0,1);
			
			double probabilitiesRatio = newProba/previousProba; 
			
			if(probabilitiesRatio>1){
				previousProba = newProba;
				numberOfSamples++;
				accept(x);
			}else{//probabilitiesRatio<1
				double threshold = gen.nextDouble();
				if(probabilitiesRatio>threshold){
					previousProba = newProba;
					numberOfSamples++;
					accept(x);
				}else{
					//reject data point
					rejected++;
				}
			}
		}while(numberOfSamples<sampleSize);
		
		FileOutputStream fos = new FileOutputStream(new File(outputFile));
		
		for(Double result: results){
			fos.write((result+"\n").getBytes());
		}
		fos.flush();
		fos.close();
		
		return rejected;
	}
	
	private void accept(double x) {
		results.add(x);
	}

	public double normal(double x, double mu, double sigma){
		double argument = (x-mu)/sigma;
		return Math.exp(-argument*argument/2)/(sigma*Math.sqrt(2*Math.PI));
	}
}
